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The Wang-Landau method is used to study the magnetic properties of the giant paramagnetic 
molecule {Mo72Fe3o} in which 30 Fe 3+ ions are coupled via antiferromagnetic exchange. The two- 
dimensional density of states g(E,M) in energy and magnetization space is calculated using a 
self-adaptive version of the Wang-Landau method. From g(E,M) the magnetization and magnetic 
susceptibility can be calculated for any temperature and external field. 
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I. INTRODUCTION 

During the past decade the field of molecular magnets 
has experienced a rapid evolutioni* 2 ^ Nowadays a vast 
variety of species can be synthesized, ranging in size from 
2 to more than 30 paramagnetic ions embedded in the 
host molecule^ The fascinating properties of these new 
materials include hysteretic behavior,— quantum tun- 
neling of the magnetization^^ magneto caloric and 
magnetostrictive effects^ 

Astonishingly, for several observables and not too low 
temperatures it is sufficient to treat magnetic molecules 
classically by applying the Heisenberg model. 12 This 
enables one to apply the powerful machinery of clas- 
sical stochastic sampling methods. The Metropolis 
algorithm 13 has become the standard tool to calculate 
statistical classical properties of these nano magnets j 14 i 15 
Recently, the Wang-Landau algorithm 16 (WL) has been 
successfully applied to various problems in statistical 
physics and biophysics both to models with discret o 16 i 17 
as well as with continuous degrees of freedom ^ 19 ! 20 ! 21 
In this context the WL exhibits a superior feature in 
comparison to the Metropolis algorithm. Since the WL 
calculates the density of states (DOS), one can estimate 
thermodynamic observables such as the free energy and 
entropy for all temperatures using just one single simu- 
lation of the density of states. 

However, the calculation of the DOS in energy space 
g(E) permits only to calculate thermodynamic properties 
as a function of temperature and at zero magnetic field, 
e.g. the zero-field specific heat, but not observables such 
as the magnetic susceptibility at non-vanishing field. 22 
For studying the magnetic properties of a system at arbi- 
trary external magnetic field one has to calculate a joined 
DOS g(E, M) in energy and magnetization space. Once 
g(E, M) is known, properties like magnetization M and 
magnetic susceptibility \ can be calculated at any tem- 
perature and any external magnetic field using again just 
one single simulation of the density of states j 21 i 23 

In this article we introduce a self-adaptive version of 
the WL which allows to calculate the two-dimensional 
DOS g(E, M) using a discrete binning scheme in the 



continuous energy and magnetization space. We demon- 
strate that the WL is capable to calculate efficiently all 
thermodynamic properties of rather large spin systems 
using the example of the magnetic Keplerate molecule 
{Mo72Fe3o}^ The dependencies of the thermodynamic 
observables both on temperature and external magnetic 
field can be obtained by only one single simulation. We 
conclude by discussing the problems arising at low tem- 
peratures and high magnetic fields. 



II. MODEL AND COMPUTATIONAL 
METHODS 

In the Keplerate molecule {Mo72Fe3o} 30 iron (III) 
ions (s = 5/2) occupy the vertices of a perfect 
icosidodecahedron^i see Fig. [TJ The spins are coupled 
with their nearest neighbors by an isotropic and antifer- 
romagnetic coupling of strength J/ks = 1.566 K. The 
spectroscopic splitting factor is g = 1.974. 12 




FIG. 1: Geometrical structure of an icosidodecahedron. For 
the molecule {Mo72Fe3o} the vertices represent spin sites and 
the edges represent nearest neighbor interactions. 

We write the Heisenberg Hamiltonian as 

h = j ]T s m -s n + 5MB 5 (z) -]Tsi z) , (1) 

(m,n) n 

whereas (m, n) directs that the sum is over distinct 
nearest-neighbor pairs, is an external magnetic field 
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in z-direction, /ig is the Bohr ma gneton, an d S denotes 
classical vectors of length |S| = y (s(s + 1))^ We use a 
self-adaptive 26 scheme of the WL to calculate the two- 
dimensional DOS g(E, M), where E denotes the Heisen- 
berg energy of a given spin configuration of Hamiltonian 
(PQ) without external magnetic field, i.e. 

E = J ^ S m • S n , (2) 

<m,n> 

and M is the magnetization in z-direction, which is de- 
fined as the sum over the z-components of all classical 
spin vectors: 

M = Y J S ( n z) - (3) 

n 

In contrast to the Metropolis algorithm, where the accep- 
tance probability of a new generated state is determined 
by min{exp[— (Ej — £^)/(/cbT)], 1}, the WL is character- 
ized by an acceptance ratio mm{g(Ej)/ g(Ei), 1}, where 
Ei and Ej refer to energies before and after the transi- 
tion. The original WL has been applied to models where 
the energy assumes only discrete values i 16 i 17 Since the 
Heisenberg model consists of classical spin vectors which 
are continuously orient able, the possible energy and mag- 
netization values are real numbers between the minimum 
and maximum energy and magnetization, respectively. 
Thus, we discretize the continuous energy and magne- 
tization range by the introduction of bins i 21 i 22 We have 
chosen bins of uniform width of ( AE) //cb = 4 K in energy 
and AM = 0.6 • |S| in magnetization range. To speed up 
the simulation we divide the total energy range in eight 
overlapping intervals and follow the recipe in Ref. [27| to 
avoid boundary effects. A priori it is not known, which 
magnetization bins are accessible by the random walker 
in a certain energy interval. A striking feature of the WL 
is, that one does not have to know anything about the 
DOS one wants to calculate. Following the original WL, 
the random walker biases itself to explore the accessible 
energy and magnetization space. To do so, we have cho- 
sen the following procedure: At the beginning we perform 
a trial run for the desired energy interval. We introduce 
an initial DOS gi n itiai(E, M) = 1 and a reference his- 
togram RH(E, M), and perform spin flips according the 
original WL procedure accepting and neglecting new spin 
configurations for a defined number of Monte- Carlo steps 
n fnitiai- One Monte-Carlo step corresponds to a single 
spin tilt event. Each time a bin (Ei,Mi) is visited the 
corresponding entries in RH(Ei, Mi) = RH(Ei, Mi) + 1 
and ginitiai(Ei,Mi) = f ■ g initial (E i, M { ) are updated, 
whereas / is the initial modification factor. Following 
this recipe the random walker triggers itself to explore 
the accessible energy and magnetization space for a pre- 
defined energy interval. 

After nf^^ ial steps have been performed we stop the 
initial run. Now we continue with the normal WL for the 
same energy interval. As an initial guess for the DOS we 
use g initial {E, M). A new state is only accepted if the 



-300 



-400 




-0.4 0.4 

M / (30- |S|) 



FIG. 2: Low-energy part of the two dimensional map for the 
reference histogram with RH(E, M) ^ for the molecule 
{Mo72Fe3o}. RH(E,M) ^ defines the accessible energy 
and magnetization space. The grid shows the bins in energy 
and magnetization space. 



corresponding entry in RH(E, M) is not zero, compare 
Fig. [21 In other words, we accept only states correspond- 
ing to bins in the energy and magnetization space which 
have already been visited during the initial run, other- 
wise the new generated state is neglected according to 
Ref. [27]. In addition, the relative flatness of the accu- 
mulated histogram is only checked after all valid entries 
in the reference histogram have been visited at least one 
hundred times. Of course the total runtime of the simula- 
tion as well as the accuracy of the finally calculated DOS 
are very sensitive to a good choice of nf n g ial . If n%g al 
is chosen very large, the random walker explores bins at 
the boundaries of the accessible energy and magnetiza- 
tion space, which are rarely visited. Thus the DOS for 
these entries is small and does not contribute significantly 
to the partition function. Nevertheless, sampling of these 
states increases the total runtime of the simulation since 
a flat histogram is to be accumulated during the simula- 
tion. On the other hand, if nf^g ial is chosen too small, 
the energy and magnetization space is sparsely explored 
and many bins of the accessible energy and magnetization 
space are not visited, resulting in an inaccurate estimate 
of the DOS. Thus, the question is to find a good tradeoff 
between runtime and accuracy of g(E,M). To estimate 
a reasonable value for nf^g ial for each energy interval, 
we have chosen different values and counted the number 
of visited distinct bins in RH(E, M) after finishing the 
initial run. It turns out that the number of visited bins 
shows some saturation behavior. After a critical num- 
ber of steps the visited energy and magnetization space 
does not grow significantly any more. Thus, after this 
critical number of steps has been reached we consider 
the accumulated RH(E, M) to be a good estimate for 
the successive WL run. As the initial modification fac- 
tor / we have chosen a rather large value of f start = e 4 
and reduce / in large steps ^ After finishing the initial 
run we perform 14 steps decreasing / according to the 
recipe fi+\ = \J~f% resulting in a final modification factor 
of 1.000000015. After g(E,M) has been obtained, the 
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magnetization M(T, B) as a function of temperature T 
and external magnetic field B can be calculated from 



(M(T, B)) 



Z E , M Mg(E,M)eM-H/k B T) 
Z EM 9(E,M)exp(-H/k B T) ' 



(4) 



The differential magnetic susceptibility x = d(M)/dB 
can be computed by using the equation 



X(T,B) 



1 



k B T 



[{M 2 (T,B)) - (M(T, B)) 2 



(5) 



While without external magnetic field it is sufficient 
to calculate only the DOS for negative energies, since 
states with positive energies due to their Boltzmann fac- 
tor exp(— E/k^T) practically do not contribute to the 
partition function, in the case of an external magnetic 
field also energy bins with a positive energy become im- 
portant, since they might be low- lying in the case of an 
applied external field. Thus it is required to calculate 
the DOS over the full energy range of the system, if 
the properties at high magnetic fields are studied. The 
exact ground state energy for {Mo72Fe3o} is known to 
be E /kB = —412.125 K? 29 which is close to the lowest 
accessible energy bin of E min /kB = —406 K. The clas- 
sical ground state of the molecule is charactereized by 
relative angles of 120° between nearest-neighbor spins. 29 
The possibility of generating such a spin configuration 
is practically zero, and thus an effective sampling of the 
DOS near the ground state energy is not possible. The 
same is valid for the highest energy states of the molecule, 
where the probability of generating a spin configuration 
with all spins pointing in the same direction is unlikely. 
The highest energy bin visited during the initial run is 
Bmax/^B = 786 K. As the criterion of flatness we used 
a relatively small value of 0.5, and to further decrease 
the statistical error we perform four runs, calculate the 
magnetic properties for each of these runs and compute 
the mean value and mean standard deviation afterwards. 
On a personal computer (Intel, Xeon, 2.60 GHz) the sam- 
pling of the complete DOS took about 120 hours of cpu 
time. 

For a comparison of accuracy we performed extensive 
Monte-Carlo simulations using the Metropolis algorithm. 
We perform 30 • 10 6 spin tilt trials to let the system reach 
equilibrium at a defined external magnetic field and tem- 
perature. Afterwards for additional 30 • 10 6 tilting trials 
the thermodynamic properties are computed. 
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FIG. 3: (Color online) Joint DOS \n[g(E,M)} for the mag- 
netic molecule {Mo72Fe3o}. 



the classical density of states grows very rapidly by or- 
ders of magnitude which explains why it is difficult to 
obtain very accurate results at the boundaries. 
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III. RESULTS AND DISCUSSION 



FIG. 4: (Color online) Magnetization M as a function of tem- 
perature T and external magnetic field B. 



The calculated joint DOS g(E, M) for the magnetic 
molecule {Mo72Fe3o} is presented in Fig. [3l One no- 
tices that the low-energy boundary E m i n (M) assumes a 
parabolic shape as observed for various systems i 3Q i 31 The 
region with dg(E, M)/dM = at high energies indicates 
the global rotational symmetry of the system. One also 
notices that in the vicinity of the phase-space boundary 



The limited accuracy at the energy and magnetization 
space boundaries looses its significance with increasing 
temperature. This is demonstrated by evaluating the 
magnetization as a function of temperature T and ex- 
ternal magnetic field B. Figure [H shows the behavior of 
the magnetization in the relevant temperature and field 
range. The simulated magnetization does not show any 



4 



significant statistical fluctuations; it compares nicely to 
the result of a Metropolis sampling. 



100 











— WL, T = 2 K 

— WL, T = 4 K 
° Metropolis, T = 2 K 
n Metropolis, T = 4 K 



10 



B[T] 



30 



40 



FIG. 5: (Color online) Comparison of the magnetization M as 
a function of external magnetic field calculated from g(E, M) 
(solid lines) with computed data using the Metropolis algo- 
rithm at T = 2 K (open circles) and 4 K (open squares), 
respectively. Statistical sampling errors are smaller than the 
used symbols and line widths. 

Figure [5] compares the magnetization at T = 2 K and 
4 K as a function of field for the WL and the Metropo- 
lis simulations. As can be inferred from the figure, at 
a temperature of the order of the exchange coupling the 
WL reaches the same accuracy as the Metropolis algo- 
rithm, but the latter has to be performed for each pair 
of variables (T,B), whereas with the WL the density of 
states has to be sampled only once in order to obtain the 
complete function (M(T, B)). 
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FIG. 7: (Color online) Comparison of the magnetic suscepti- 
bility x as a function of external magnetic field calculated 
from g(E,M) (solid lines) with computed data using the 
Metropolis algorithm at T — 2 K (open circles) and 4 K 
(open squares), respectively. 



mum at about £> sat /3, 14 compare Fig.0 Clearly one can 
see that the dip in suszeptibility at about 5^/3 van- 
ishes with increasing temperature^ The spurious peak 
at T = and B = reflects the difficulty to obtain the 
density of states in close vicinity of the ground state. 
A comparison with Metropolis simulations is shown is 
Fig. The results obtained using both methods agree 
well. The fluctuations in the WL data are due to the 
missing bins in g(E, M) at the energy and magnetiza- 
tion space boundaries, since these features are apparent 
in all four WL runs. The statistical sampling errors are 
given by the size of the symbols and by the line widths 
used in Fig. [7| 




T[K] 



B[T] 



FIG. 6: (Color online) Magnetic susceptibility x as a function 
of temperature T and external magnetic field B. 




FIG. 8: (Color online) Specific heat C as a function of tem- 
perature T and external magnetic field B. 



The differential susceptibility, Figs. [6] and [7J is a sec- 
ond derivative of the partition function. Therefore, it will 
magnify inaccuracies of the simulated density of states. 
Figure [6] displays the susceptibility as a function of T 
and B. It can be seen that for temperatures higher than 
J/ks the behavior is rather smooth, but for lower tem- 
peratures a spiky behavior is observed which results from 
statistical fluctuations at the energy and magnetization 
space boundary. At T « J/ks these fluctuations are still 
visible, but on average much smaller than the real mini- 



The last function we want to discuss in this paper is 
the specific heat C as a function of temperature T and 
external magnetic field B. As can be seen in Fig. [8] the 
specific heat does not show too strong fluctuations at low 
temperature and fields below saturation. It seems that 
this observable is more robust against statistical fluctu- 
ations of the density of states than the susceptibility. In 
order to complete the discussion we display the specific 
heat also for magnetic fields above the saturation field 
where it shows strong fluctuations. They reflect magni- 
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fled inaccuracies at the high energy boundary of g(E, M) 
and thus are spurious. 

Summarizing, one can say that the proposed self- 
adaptive version of the Wang-Landau algorithm can ef- 
ficiently generate the density of states g(E, M) of a 
rather large spin system, such as the magnetic molecule 
{Mo72Fe3o}. The obvious advantage is that with one sim- 
ulation of g(E, M) many thermodynamic observables can 
be evaluated as function of temperature and applied mag- 
netic field. Statistical fluctuations are apparent in the 
second derivatives of the two-dimensional DOS, namely 
the specific heat and magnetic suszeptibility at low tem- 
peratures. These fluctuations can be minimized at the 
cost of computational time, by either increasing the num- 
ber of bins or by a more strict flatness criterion. Recently, 
it was demonstrated that a fitting of the one-dimensional 
DOS for a sampled energy interval by higher order poly- 
nomials results in a less fluctuating specific heati^ Nev- 



ertheless, it is stated, that this fitting procedure fails if 
one tries to extrapolate the DOS outside the sampled 
energy interval using the fitting function. Thus, further 
improvements are still needed to provide an effective sam- 
pling at the boundaries of g(E, M). Due to the fact that 
the density of states varies by many orders of magnitude 
it is clear that close to the boundaries of g(E, M) the 
statistics becomes poorer. 
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